%%%% This code is to produce table 6 panel A %%%%%%%%%%%%
clear
cd ('F:\home\Dropbox\Mortgage_company\JobMarketProject\Replication_Package\financing_competitors_upload\code\');

%% Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 41);
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["Year", "MSA", "LoanPurpose", "rho0", "rho1", "msa_tot_raw", "income", "pop", "income_per_cap", "mean_income", "sd_income", "rate_mb", "rate_wb", "rate_sb0", "rate_sb1", "County", "State", "CensusTract", "N_sb0", "N_sb1", "N_mb", "N_wb", "sb0_count", "sb1_count", "mb_count", "wb_count", "msa_tot_count", "outside", "scaler", "market_size", "sb0_share", "sb0_ind_share", "sb1_share", "sb1_ind_share", "mb_share", "mb_ind_share", "wb_share", "wb_ind_share", "outside_share", "sbw", "sow"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
input = readtable("..\data\input.csv", opts);
clear opts
%% parameters 
F= 18380*354500; % MSA* average loan size
alpha = 1.14;
eta = 0.09;
k = 3.21; %kappa 
qb = 2.43; 
qn = 1.67;
sigma = 0.30; % 0.50 from cross-section
%%
psiw_output = [];
psim_output = [];
for i = 1:size(input,1)
    i
    diff = 1000;
    diff2 = 1000;
    iter = 0;
    delta = 0.1;
    
    rhob = input.rho1(i);
    rhoo = input.rho0(i);
    rb = input.rate_wb(i);
    sb = input.wb_ind_share(i); % bank mortgage market share 
    Nn = input.N_sb0(i)+input.N_sb1(i);
    Nb = input.N_wb(i);
    Nm = input.N_mb(i);
    sbw = input.sbw(i);
    sow = input.sow(i);
    snb = input.sb1_ind_share(i);
    
    psiw_initial=3;
    while diff>10^(-9) && iter<10000000
        iter = iter+1;
        
        rhobl = rhob+delta;
        % calibration: psiw; find drb dsb dsnb first
        psim_initial = rb-1/(alpha*(1-sb))-eta*(rhob-psiw_initial)*(Nn*sbw*snb)/(1-sb); 


        [~, ~, ~,snbl1,~,~,~,~,rbl1,sbl1,~,~] = Solve_profit_wbank_other(psim_initial,psiw_initial,alpha,k,Nb,Nn,Nm,qn,qb,eta,rhoo,rhob,rhob,sow,sbw);
        [~, ~, ~,snbl2,~,~,~,~,rbl2,sbl2,~,~] = Solve_profit_wbank_other(psim_initial,psiw_initial,alpha,k,Nb,Nn,Nm,qn,qb,eta,rhoo,rhob,rhobl,sow,sbw);
               
        drb = (rbl2-rbl1)/delta;
        dsb = (sbl2-sbl1)/delta;
        dsnb = (snbl2-snbl1)/delta;
        psiw = rhob-(snbl1/(-dsnb*(1+sigma*(1-sbw)/(snbl1*(1-snbl1))))+(drb*sbl1+(rbl1-psim_initial)*dsb)*(eta*Nn*sbw*(-dsnb*(1+sigma*(1-sbw)/(snbl1*(1-snbl1)))))^(-1));

        diff = abs(psiw-psiw_initial);
        psiw_initial = psiw;
        %psiw_initial = rhob-((rb-1/(alpha*(1-sb))-psim)*(1-sb)/(eta*(Nn*sbw*snb)));
        psim_initial = rb-1/(alpha*(1-sb))-eta*(rhob-psiw_initial)*(Nn*sbw*snb)/(1-sb); 

    end
    psiw_output(i) = psiw_initial;    
    psim_output(i) = psim_initial;

end 

       
psiw = mean(psiw_output)
psiw_str = std(psiw_output)/sqrt(size(psiw_output,2))
psim = mean(psim_output)
psim_str = std(psim_output)/sqrt(size(psim_output,2))
%% data 
rhob_data =  mean(input.rho1)
rhoo_data = mean(input.rho0)
rb_data = mean(input.rate_wb)
rm_data = mean(rmmissing(input.rate_mb))
rnb_data = mean(rmmissing(input.rate_sb1))
rno_data = mean(rmmissing(input.rate_sb0))
sb_data = mean(input.wb_ind_share.*input.N_wb) % bank mortgage market share 
sbm_data = mean(rmmissing(input.mb_ind_share.*input.N_mb)) % bank mortgage market share 
sbw_data = mean(input.sbw.*input.N_wb)
sow_data = mean(input.sow)
snb_data = mean(input.sb1_ind_share.*input.N_sb1)
sno_data = mean(rmmissing(input.sb0_ind_share.*input.N_sb0))
outside_data = mean(rmmissing(input.outside_share))
%% Model prediction 

rhoo = mean(input.rho0);
Nn = mean(input.N_sb0+input.N_sb1);
Nb = mean(input.N_wb);
Nm = mean(input.N_mb);
[rno,sno,rnb,snb,rm,sbm,rb,sb, rhob, sbw, sow,~,~] = Solve_local(F,sigma,psim,psiw,alpha,k,Nb,Nn,Nm,qn,qb,eta,rhoo);
rhob
sbw*Nb
sow 
rb
rm 
rnb
rno
sb*Nb
sbm*Nm
snb*Nn*sbw*Nb
sno*Nn*sow